STABILISED FINITE ELEMENT METHODS FOR 

NON-SYMMETRIC, NON-COERCIVE AND ILL-POSED PROBLEMS. 

PART I: ELLIPTIC EQUATIONS 

ERIK BURMAN* 

Abstract. In this paper we propose a new method to stabiUse non-symmetric indefinite prob- 
f"C^ , lems. The idea is to solve a forward and an adjoint problem simultaneously using a suitable stabilised 

finite element method. Both stabilisation of the element residual and jumps of certain derivatives of 
f**^ ' the discrete solution over element faces may be used. Under the assumption of well posedness of the 

^vj ' partial differential equation and its associated adjoint problem we prove optimal error estimates in 

H^ and L^ norms in an abstract framework. Some examples of problems that are neither symmetric 
nor coercive, but that enter the abstract framework are given. First we treat indefinite convection- 
^ , diffusion equations, with non-solenoidal transport velocity and either pure Dirichlet conditions or 

■^^ ' pure Neumann conditions and then a Cauchy problem for the Helmholtz operator. Some numerical 

_ ' illustrations are given. 
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1. Introduction. The computation of indefinite elliptic problems often involves 
certain conditions on the mesh size h for the system to be well-posed and for the 
derivation of error estimates. The first results on this problem are due to Schatz [53]. 
(-H ' Such conditions on the mesh parameter can be avoided if a stabilised finite element 

method is used. Such methods have been propsed by Bramble et al. HO] or more 
recently the conitnuous interior penalty (GIF) method for the Helmholtz equation 
suggested by Wu et al. [351 HI]- The method proposed herein has some common 
features with both these methods, but appears to have a wider field of applicabil- 
ity. We may treat not only symmetric indefinite problems such as the (real valued) 
^ ■ Helmholtz equation, but also nonsymmetric indefinite problems such as convection- 

'^ I diffusion problems with non-solenoidal convection velocity or the Cauchy problem. 

The latter problem is known to be ill-posed in general [2] and will, in the ill-posed 
rsS I case, mainly be explored numerically herein. For all these cases we show that if the 

■ ' ■ primal and adjoint problems admit a unique solution with sufficient smoothness the 

^— s I proposed algorithm converges with optimal order. The case of hyperbolic problems is 

fT^ ■ treated in the companion paper (^. 

The idea of this work is to assume ill-posedness of the discrete form of the PDE 
and regularise it in the form of an optimisation problem under constraints. Indeed 
we seek to minimise the size of the stabilisation operator under the constraint of the 

discrete variational form. The regularisation terms are then chosen from well known 
s_j , ....... 

jrt , stabilised methods respecting certain design criteria given in an abstract analysis. 

This leads to an extended method where simultanuously both a primal and a dual 

problem are solved, but where the exact solution of the dual problem is always trivial. 

The aim is to obtain a method where possible discrete non-uniqueness is alleviated by 

discrete regularisation with a non-consistency that can be controlled so that optimal 

convergence for smooth solutions is obtained. The method is also a good candidate for 

cases where the solution to the continuous problem is non- unique, but that is beyond 

the scope of the present paper. 

In spite of the lack of coercivity for the physical problem, the discrete problem 

has partial coercivity on the stabilisation operator. A consequence of this is that, 

depending on the kernel of the stabilisation operator a unique discrete solution may 
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often be shown to exist independently of the underlying partial differential equation. 
This can be helpful when exploring ill-posed problems numerically or when measure- 
ment errors in the data, may lead to an ill-posed problem, although the true problem 
is well posed. 

An outline of the paper is as follows, in Section[2]we propose an abstract method 
and prove that the method will have optimal convergence under certain assumptions 
on the bilinear form. Then in Section [3] we discuss stabilised methods that satisfy 
the assumptions of the abstract theory with particular focus on the Galerkin least 
squares method (GLS) and the continuous interior penalty method (CIP). Three 
example of applications are given in Section |4l two different non-coercive transport 
problems in compressible flow fields and one elliptic Cauchy problem. Finally in 
Section [5] the accuracy and robustness of the proposed method is shown by some 
computations of solutions to the problems discussed in Section |4l In particular we 
study the performance of the approach for some different Cauchy problems of varying 
difficulty. 

2. Abstract formulation. Let fi be a polygonal/polyhedral subset of M'^. For 
simplicity we will reduce the scope to second order elliptic problems, but the method- 
ology can readily be extended to indefinite elliptic problems of any order, providing 
the operator has a smoothing property. 

We let V,W denote two subspaces of H^{fl). The abstract weak formulation of 
the continuous problem takes the form: find u Cz V such that 

a{u,v) ^ {f,v), '^veW (2.1) 

with formal adjoint: find z ^W such that 

a(w, z) = (g, w), Vw e V. (2.2) 

The bilinear form a(-, •) : V x 14^ — !■ R is assumed to be elliptic, but neither symmetric 
nor coercive. We denote the forward problem on strong form Cu = f and the adjoint 
problem on strong form C* z = g. Suitable boundary conditions are integrated either 
in the spaces V, W or in the linear form. 

We assume that both these problems are well posed and that the geometry and 
data are such that the smoothing property holds 

|w|h2(o) < Ca,fi||/||, \z\H'^(n)<Ca,n\\9\\- (2.3) 

We will frequently use the notation a ^ b for a < Cb with C a constant depending 
only on the mesh geometry and physical parameters of order one (or irrelevant in the 
context). We will also use a ^ b for a <b and b < a. Indexed constants Cxy will have 
similar dependence on the variables xy at each occurence but can differ by a constant 
factor. 

The L^-scalar product over some X C M'' is denoted {■,-)x and the associated 
norm || • \\x, the subscript is dropped whenever X — fl. We will also use (•, •}y to 
denote the L^-scalar product over Y C R'^"^ and {■,-)h the element- wise L^-norm 
with the associated broken norm || ■ \\h- 

2.1. Finite element discretisation. Let {Th}h denote a family of quasi uni- 
form, shape regular triangulations Th '■— {K}, indexed by the maximum triangle 
radius h. The set of faces of the triangulation will be denoted by J- and Tint denotes 

2 



the subset of interior faces. Now let X!^ denote the finite element space of continuous, 
piecewise polynomial fuctions on Th, 

X,f := {Vh e H\n) : Vh\K e Pk{K), ^K G %}■ 

Here Pfc (K) denotes the space of polynomials of degree less than or equal to fc on a 
triangle K. 

We let TTi denote the standard L^-projection onto XJl and ih '■ C"(r2) t-^ XJl the 
standard Lagrange interpolant. Recall that for any function u G (F U W) Ci H''~^^{Q) 
there holds 

\\u - ihu\\ + h\\V{u - ihu)\\ + h^\\D^{u - ihu)\\h < c,h''+^\u\H>^+^n)- (2.4) 

and under our assumptions on the mesh, similarly for tt^. We propose the following 
finite element method for the approximation of (|2.ip . find {uh,Zfi) G 14 x Wh such 
that 

(2.5) 
ah{wh,Zh) - Sp{uh,Wh) = -Sp{u,Wh), 

for all {vh,Wh) € Wh x Vh- The spaces Vh and Wh are typically both chosen as 
X^, possibly with some additional constraint. If both the primal and the adjoint 
variable are approximated in the same spaces, we use the notation Vh for both. The 
bilinear form ah{-, •) is a discrete realisation of a(-, •), typically modified to account for 
the effect of boundary conditions, since in general Vh ^ V and Vh ^ W. The penalty 
operators Sa{-,-) a-nd Sp{-,-) are symmetric stabiHsation operators and associated with 
the adjoint and the primal equation respectively. The rationale of the formulation 
may be explained in an optimisation framework. Assume that we want to solve the 
problem, find Uh G Vh such that 

ah{uh,Vh) ^ if.Vh), yvheVh, 

but that the system matrix corresponding to ah{uh,Vh) has zero eigenvalues. The 
discrete system is ill-posed. This often reflects some poor stability properties of the 
underlying continuous problem. The idea is to introduce a silection criterion for the 
solution, in order to ensure discrete uniqueness, measured by some operator Sp{uh, Vh)- 
This can include both stabilisation (regularisation) terms and the fitting of the com- 
puted solution to measurements or some goal function. Interestingly, as we shall 
see, finite element stabilisation terms are the only regularisation that is needed. The 
formulation then writes, find Uh, z/i G V/i x Wh stationary point of the Lagrangian 



-Sp{Uh ~U,Uh-U) - -i 



L{uh, Zh) := 7;Sp{uh -u,Uh-u) - -Sa{zh, Zh) - ah{uh, Zh) + if, Zh). (2.6) 



The saddle point structure of the Lagrangian has been strengthened by the addition 
of the regularizing term —\sa{zh, Zh). We may readily verify that 



dL 

duh 

and 

dL 



■^;—{wh) = Sp{uh - u,Wh) - ah{wh,Zh) 



r. (vh) = -ah{uh,Vh) - Sa{zh,Vh) + if,Vh). 



It follows that (12.51) corresponds to the optiniality conditions of 

Observe that the second equation of (I2.5P is a finite element discretisation of the 
dual problem (|2.2[) with data 5 = 0. Hence the solution to approximate is z = 0, how- 
ever Zfi will most likely not be zero, since it is perturbed by the stabilisation operator 
acting on the solution Uh, which in general does not coincide with the stabilisation 
acting on u. 

We will assume that the following strong consistency properties hold. If u is the 
solution of (|2.ip then 

ahiu,f)^{Cu,ip) ioT al\ipeW + Wh (2.7) 

and if z is the solution of (|2.2I) then 

a,,((/), z) = ((/), C*z) for all G F + V^. (2.8) 

As a consequence the following Galerkin orthogonalities hold 

ah[u " Uh.Vh) = Sa{zh,Vh) and ah{wh,z - Zh) = Sp[u-Uh,Wh). (2.9) 

The bilinear forms Sa(-, ■), Sp(-, •) are symmetric, positive semi-definite, weakly consis- 
tent, stabilisation operators. For simplicity we will always assume that u is sufficiently 
regular so that Sp{u,Vh) is well defined, i.e. the stabilisation is strongly consistent. 
The analysis usin weak consistency of the stabilisation only is a straightforward mod- 
ification. The semi-norms on Vt and Wh associated to the stabilisation is defined 

by 

\xh\Sy ■= Sy[xh,Xh)^, y^a,p. 

The rationale of this formulation is that the following partial coercivity is obtained 
by taking Vh = Zh and Wh = u^, 

\zh\l^ + \uh\l^ = {f,Zh) -Sp{u,Uh). (2.10) 

We assume that there are interpolation operators ny ■ V —i' Vh and ttw ■ W — >■ 
Wfi, satisfying (|2.4[) and also the following continuity relations for all v,w,y £ H^{il) 
and for all Xh £ Wh 

ah{v -TTvV.Xh) < \\v ~ Tlvv\\j^{Ca\Xh\Sa. +^{h)\\xh\\) (2-11) 

and 

ah{u-Uh,y-Trwy) < \\y-'rrwy\\*{ca\\u-Trvu\\c+Ca\u-Trvu\sp+e{h)\\u-Uh\\) (2.12) 

where || • ||-i-, || • ||* and || ■ ||£ are semi- norms to be defined, satisfying the approximation 
estimates 

\\v-TTvv\\c+\\v-TTvv\\ + + \v-Trvv\s^ < Ca,jh''\v\H>^+im), Vw e Vr\H''+\n), (2.13) 



Iw-TTwwW^ + lw -Trww\s^ <Ca,^h''\w\Hk+im) \fw eW r]H''+^{n), (2.14) 



i(n), 



and 



It^wwIs^ < Ca,ih\w\H2{n), Vw e wnH'^{Q), and \nvv\sj, < Ca,7/i|'y|i/2(o) Vv G VnH'^{fl)^ 

(2.15) 
where Ca^-y depends on the form a(-, •) and a stabilisation parameter 7. 
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2.2. Convergence analysis for the abstract method. We first prove that, 
assuming that the discrete solution (u^, z/j) exists, the stabihsation semi-norm of the 
discrete error is bounded by one term that converges to zero at an optimal rate and 
another non-essential perturbation that can be made small. 

Lemma 2.1. Assume that the solution of (|2.5p exists, that the solution of (|2.ip is 
smooth, that the forms of ()2.5|) and the operators ixy , ttw are such that (|2.9|) . (|2.1ip 
and ()2.13|) are satisfied. Then there holds 

\Trvu-Uh\sp + \t^wz - Zh\sa < CQ,^,<:/i''|w|^fc+i(n) + e{h)\\zh\\, 

1 

where Ca,-y,e '■= '^a.-yi'^ + c'^)'^ , with Ca and c a, -J defined by (I2.11|) and p.l3|) respectively. 
Similarly, if Sp{u,Wh) ~ 0, there holds 

\Uh\Sp + \zh\s^ < (Ca,7 + Ca,^,e)h''\u\Hk+l(^n^ + e{h)\\zh\\. 



Proof. For the first inequality let ^h = ttvu — Uh and Ch = ttwz — Zh- By the 
definition (12.51) there holds 



l^ftlSp + ICftlL = Sp{S,h,^h)+Sa{C,hXh) ^ah{^h,Ch)+SaiCh,Ch)-ah{^h,Ch) + Sp{^h,^h)- 

Using now the Galerkin orthogonality of a;i(-, •), p.9p . we have 

IChllp + IC/i||„ = a/i(7ryu-u, 00 + SaiTTwzXh) -ahi^h,'!rwz- z) + Sp{-Kvu - u,^h)- 
Observing that z — t:\yz — this reduces to 

IC/ilsj, + IC/J|„ = ai^{nvu~u,Ch) + Spijiyu ~ u,S,h)- 
We conclude by applying the continuity (|2.1ip 

If/ilsp + IOi||„ < \\u - TTvu\\ + {ca\C,h\s^ +e{h)\\zh\\) + \u - ■nvu\s^\^h\s^ 
followed by 

laiL + laiL < {cl + l)\\u - -.vuWl + eihfWz^f + \u- ^vu\l^ 

.,7(1 + (^a)^ l"l//'« + i(0) 



<cl^{l+cl)h^%\l,^,,^.+e{hf\\zuf. 



The second result follows by adding and subtracting Tryw, observing that here -k^z = 
0, applying a triangle inequality and then (|2.13p on \'n:vu\sp — \''^vu — u\sj,- D 
We may now prove the main result which is optimal convergence in the L^ and the 
H^ norms. 

Theorem 2.2. Assume that (|2.ip and (|2.2p are well-posed with exact solutions 
u and z satisfying (j2.3p . Assume that the forms of (j2.5p and the operators Try, tt\y 
are such that (j2.9p - (|2.15p are satisfied and that h is so small that 

Ca,j.Qhe{h)<-, (2.16) 

6 

wif/i Ca.-y^n depending mainly on the constants of the inequalities (|2.3p and ()2.13p - 
p.lSp . T/ien f/ie discrete solution Uh, Zh exists and satsifies 

\\u~Uh\\ +h\\V{u-Uh)\\ + \\zh\\ < Ca,n,'yh''^^\u\Hk+ii^n) 



and in particular 

\\u - UhW + h\\Viu - Uh)\\ + \\zh\\ < Ca,nnh^\\f\\. (2.17) 

Proof. Let (p be the solution of (|2.2p with g — u — ut and ip the solution of (|2.ip 
with f = Zfi- By (12.31) there holds 



WvWH-^in) < Ca,n\\u - Uh\\ and ||V'll_ff2(o) < Ca,n\\zh\\- 

By definition of the primal and dual problems and by ((2?7t . ((2?8l) . ((2^ . (|2.1ip and 
(I^TT^ there holds, 

\\u - M/ilP + ||zh|p = (m - Uh,C*ip) + {Cijj,Zh) = a;i(w - u;i,(y5) +ah(V', 2/1) 

= ah{u-Uh,lf-TTwf) +Sa{Zh,TTwV') 

+ ahiijj - TTvi', Zh) - Sp{u - Uh^TTvi') 
< (Ca\\u - TTyuWc + Ca\Uh -TTyUJSp + e{h)\\u - Uh\\)\\ip - TTiyvW* 

+ ica\zh\s^ + e{h)\\zh\\)U - nv^U 
+ IzhlsjT^w^lsa + \uh - "IspkvV'Isp- 

First we observe that by ^J^, (pji)) and (P3| 

eih)\\u- Uh\\\\ip - TTw^W* + eih)\\zh\\\\ip -7ryVll + 

< Ca,j.nh{e{h)\\u - M,,f + e(/i)||z,,||2). 

Then by Lemma 12.11 and approximation we have 

IzhlsjTTWfils^ + \uh - w|5pkv0|Sp 

< ((Cq,7 + Ca,.y.,e)h''\u\Hk + l(^n) + fX^^)\\zh\\)Ca,-y .nHWu - Uh\\ + \\zh\\). 

Using the two previous bounds and a arithmetic-geometric inequality we have 

(1 - 3ca,j,nhe{h))i\\u - w,,f + Wzhf) 

< Ca,'fh''^^\u\Hk+i(^n){\v\H'^in) + 1-01^2(0)). 

Using equation (|2.3p . the result for the L^-norm follows provided h satisfies (|2.16p . 
The result for the iJ-'^-norm follows using a global inverse inequality on the discrete 
error and then the L^-norm error estimate. 

||V(u - Uh)\\ < ||V(w - 7ryM)|| + \\\7{ttvu ~ Uh)\\ < h''\u\H<^+i(n) + /i"Nkyu - Uh\\. 

The existence of a unique solution to (|2.5p is a consequence of (I2.17P . Well-posedness 
of (|2.1I) means that f — implies u = 0, but then by (|2.17p Uh = Zh = 0, which shows 
that the matrix is invertible. D 
The optimal convergence of the stabilisation terms follows. 

Corollary 2.3. Under the assumptions of Lemma \2.1\ and Theorem \2.2\ there 
holds 

{ttvu - Uh\sp + Vwz - z/ilsa ^ Cs,Ji^\u\iik^i(^-^ + 0{h^^^). 



Proof. Immediate consequence of Lemma 12.11 and Theorem 12.21 D 

Remark 1. The need to control a low order contribution of the dual solution 

usually comes from oscillation of data. Either in the form of stabilisation terms that 

do not account for oscillation within the element or error in the numerical quadrature. 

In case a Gardings inequality holds for (j2.5p and Sa{-, •) = Sp{-, •) the i/^-error can be 

recovered without using inverse inequahties as stated below. 

Corollary 2.4. Assume that for the bilinear form a(-, •) there exists A € M such 

that 

llVuhiP - A||d/J^ < ah{vh,Vh) + Sp{vh,Vh) 
and that Sa{-, •) = Sp{-, •). Then 

l|V(u-Uh)|| </i1u|ff.+i(a). 



Proof. Similar to proof of Lemma 12.11 and therefore only sketched. Let S^h '■= 
■Kyu — Ufi- It follows by Gardings inequality that 

iivaf < a/.(a,a) + AiiaiP + 5p(a,a)- 

Using Galerkin orthogonality we have 

ah{£.h,£,h) = ah{Trvu-u,^h) + Sa{zh,^h) 

and the rest follows as in Lemma [2.11 bv (|2.11l) . (|2.13p and using the known conver- 
gences of Lemma 12.11 and Theorem 12.21 D 

3. Stabilisation methods. To fix the ideas let £ be a second order elliptic 
operator on conservation form, 

Cu:^ -^Au + V ■{l3u) + cu. (3.1) 

Here fj, G M"*", (3 e [C^(Jl)]^ is a non-solenoidal velocity vectorfiel and c G (7^(17). A 
necessary condition for the continuities (|2.1ip and (|2.12p are that the the stabilisation 
operators satisfy 

inf y^\\h{jCuh-Wh)\\K + \\hitiHVuh-nFJ\\jr. <Sp{uh-u,Uh-u) (3.2) 



inf y^\\hiC*Zh~Wh)\\l + \\h-2^l^Vzh■nFl\\%,^^<Saizh,Zh) (3.3) 

at least up to a non-essential low order perturbation. Note that taking Wh ^ f in 
the first term in the left hand side results in a least squares term on the residual over 
the element. It follows that the stabilisation relies on two mechanisms: L^-control 
of the element residual and L^-control of the gradient jump over element edges. If 
higher order differential equations are considered, jumps of higher derivatives must be 
added. The design criterion p.2p - p.3|) makes it straightforward to adapt the analysis 
below to a range of stabilisation methods, such as Galerkin least squares, orthogonal 
subscales, continuous interior penalty or discontinuous Galerkin methods. In all cases 
however the jumps of the gradient must be penalised, or an equivalent stabilisation 
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operator introduced. It therefore seems natural to consider two stabilisations in more 
detail, first the GLS-stabilisation combined with gradient penalty and then a CIP- 
stabilisation purely based on penalty on jumps of derivatives of the approximate 
solution. We introduce the stabilisation operators 



and 



where 



Sp{uh,Vh) := Sp.GLsiuh, Vh) + Scip{uh, Vh) (3.4) 

Sa{zh, Wh) ■■= Sa.,GLs{zh,Wh) + Sctp{uh, Vh) (3.5) 



s 



■p,GLs{uh,Vh) := {iGLsh Cuh,Cvh)h 



Sa,GLs{zh,Wh) := {-fGLsh"^ C.* Zh: C*Wh)h 



and 



Scip{uh,vh) -^ J2 (hFliA^^hi ■ lyvhl + h%-f2,Fl^UhllAvhi) dx. (3.6) 

Here |Vm;i||f and |AM/i||i? denotes the jump of the gradient and the Laplacian over 
the face F. Note that for smooth u, Sp{u,Vh) = Sp^GLs{u,Vh) = {LlGLsh'^Lvhjh 
showing the Petrov-Galerkin character of the GLS-method. The abstract analy- 
sis typically holds for the parameter choices 7gl5 > 0, 71, f > 0, 72, f = or 
IGLS = 0, 71, F > 0, 72, F > 0. Note that the matrix stencil for finite element methods 
remains the same for both approaches, and therefore the CIP-method seems more ap- 
pealing in this context. Eliminating the Galerkin least squares term also reduce the 
computational effort since the same stabilisation is used for the primal and adjoint so- 
lution. If on the other hand a C^-continuous approximation space is used, the jumps 
of the gradients may be omitted and the GLS stabilisation might prove competitive, 
since integrations on the faces may then be avoided. 

3.1. Galerkin- least-squares stabilisation. The GLS-method is one of the 
most popular stabilised methods. To fix the ideas we will assume that the prob- 
lems p.ip and (12. 2p are subject to homogeneous Dirichlet conditions and well-posed, 
with / S L'^{^). For the readers convenience we detail the Lagrangian p.6p in this 
particular case 

lj{uh,Zh) := -\\T^{Cuh - f)\\l + ■^Sctp{uh,Uh) 

- -\\T^C*Zh\\l - -Scip{zh,Zh) -a{uh,Zh) + (f,Zh). (3.7) 

The optimality conditions then write, find Uh, z^ G Vh H -ffp (fJ) 

aiuh,vh) + Sa{zh,Vh) = if,Vh) (3.8) 

a{wh,Zh) - Sp[uh,Wh) = -Sp{u, Wh) = -(/, TCwh)h 

for all Vh,Wh e Vh n Hg{i}). Here jGLsh'^ := t > and 71, f ^ M) 72, f = 0. We 
assume that the physical parameters are all order unity for simplicity. Observe the 



nonstandard structure of the stabilisation terms and that the formulation is consistent 
for u the exact solution of (|2.ip and z — 0. We will now prove that the assumptions 
of Proposition 12.11 and Theorem l2.2l are satisfied for the formulation (|3.8p . 
We define the following semi-norms 

\\vU:^\\v\U:^\\r-iv\\+( Y. \\fiih~h\\%] , (3.9) 



\\v\\c:^\\r--livl,^\\h-^pi-^l^u-nF\\\^^^, 

and 

kisp := ||T5/:a;||/i + Scip(a;, a;)2 and \x\s^ := ||r2/:*x||,, + Scip(a;,a;) ^ , 

defined for x £ H^{n,) + Vh- Let Try and ttw be defined by the Lagrange interpolator 
ih (or any other iJ^-stable interpolation operator that satisfies boundary conditions), 
and note that by (12.41) we readily deduce the following approximation results for 
smooth enough functions u 

\\u-TTvu\\++\\u-'Ilvu\\c + \u-TTvu\Sp + \\u-TTwu\\* + \u-1Twu\s„, < Ca,^h'' \u\Hk + l(^n) 

and, by i?^-stability of the interpolation operator, 

kvwisp < c^Mh\\v\\m(n), \T^Ww\sa < Ci,ah\\w\\H2(n), Vv,w £ H^{n). 

This shows that (|2.13p and p.l4p holds. It then only remains to show the continuities 
()2.1ip and ()2.12p . First we show the inequality (|2.1ip For the second order elliptic 
problem we note that after an integration by parts and Cauchy-Schwarz inequality, 

ahiv -TTvv.Xh) = ^ {u - TTyu, IfiVxh ■ npDp + y^(" ~ t^vu, C*Xh)K 

F<^Tm k 

< \\u- Trvu\\ + \xh\sa- 

Similarly, to prove (j2.12p we integrate by parts in the opposite direction in the second 
order operator and obtain 

ah{u -Uh,y- TTwy) = ^ {C{u - Uh),y - '^wy)K 
KeTh 

\\y - TTWyhiWu - T^VUWC + \uh - TTytilSp) 

Remark 2. Note that for the GLS-method e(h) = in ((2lT|) and ((27T2)) . This 
follows from the fact that the whole residual is considered in the stabilisation term. 
This nice feature however only holds under exact quadrature. When the integrals are 
approximated, the quadrature error once again gives rise to oscillation terms from 
data that introduces a non-zero contribution to e{h). 
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3.2. Continuous interior penalty. Since in this case we must account for 
possible oscillation of the physical coefficients we postpone the detailed analysis to 
the examples below and here only discuss the general principle. In this case we use 
IGLS = 0, 7i,F > 0, 72, F > in the general expressions for the stabihsation (I3.4|) and 
(|3.5p . The parameters 7i.F, i — 1,2 are stabilisation coefficients, the form of which 
will be problem specific and will be given for each problem below. The key observation 
is that the following discrete approximation result holds for suitably chosen ji^p in 
sc»p(-,-) (see [ZIIH]) 



\\h^ph ■ VU,, - losPh ■ Vw^f + ^ \\h{Auh - IosAuh)f < Sc^piUh,Uh). (3.10) 

K 



Here (3h is some piecewise affine interpolant of the velocity vector field /3 and los is the 
quasi-interpolation operator defined in each node of the mesh as a straight average of 
the function values from triangles sharing that node (see [S]). For example, 



E 

{K-.XiGK} 



{IosAuh){Xi) = iV^ ' 2^ Auh{Xi)\K, 



with Ni := cardjilT : Xi E K}. Using (|3.10p one may prove that 



inf \\Cuh-Vh\\h<Scipiuh,Uh)'^+e{h)\\uh\\. (3-11) 

Vh&Vh 



It immediately follows that p.2p and p.3p are satisfied. We will leave the discussion 
of (|2.1ip - (|2.14p and p. lip to the applications below, giving the explicit form for 
e(/i) for each case. Here we instead proceed with an abstract analysis, ignoring the 
contribution from weakly imposed boundary conditions and assuming that all physical 
parameters are of order 0(1). Here we choose Try and ttw as the L^-projection in 
order to exploit orthogonality to "filter" the element residual. We let || • ||+ and || • 1|* 
have the same definition as in the GLS case and define 



\u\\c := \\hCu\\h + Wh^fi^Vu ■ nFllb.„, + e(Mhll- (3-12) 



Then we proceed similarly as for GLS, but we use the orthogonality of the L^- 
projection, ignoring here the contribution from boundary terms. It then follows using 
the orthogonality of the projection that formally 



ah{v -Trvv,Xh) ^ ^ {u - 7rvu,lfi'Vxh ■ npl) p + '^{u - TTyu, C*Xh - Wh)K 

< \\u - Trvu\\ + {\xh\sa + e(^)l|a;/i||)- 
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Similarly, to prove (j2.12p we integrate by parts in the opposite direction in the second 
order operator and use the L^-orthogonality to obtain 

ah{u -uh,y- TTwy) = ^ m-u -~uh),y- '^wy)K 

+ X! {l^^^Uh■nF\,y-■Kwy)F 

= ^ {C{U - TTyu) + CiTTvU - Uh) -Wh,y ~ TTwy)K 

Ken 
+ X! {lt^^Uh-nFi,y-TTwy)p 

< \\y - Trwy\\*{\\hC{u - 7ryu)||^ + {ttvu - Uh\sp + e{h)\\TTvu - Uh\\) 

< \\y - T^WyhiWu - TTVU\\C + \Uh - TTyuISp + e{h)\\u ~ Uh\\) 

The last inequality follows by adding and subtracting u in the last norm in the right 
hand side to obtain e(h)\\TTvu — u + u — Uh\\- This term is then split using a triangular 
inequality and the approximation error integrated in the || • ||£ term. To use the 
L^-projection in this fashion we must impose the boundary conditions weakly so that 
the boundary degrees of freedom are included in Vh- For the GLS- method one has 
the choice between weak and strong imposition of boundary condition. In the next 
section we will discuss how weakly imposed boundary conditions are included in the 
formulation. 

3.3. Imposition of boundary conditions. The imposition of weak boundary 
conditions in this framework uses Nitsche type formulations. However they differ from 
the standard Nitsche boundary conditions in several ways: 

- Both Dirichlet and Neumann conditions are imposed using penalty. 

- There is no lower bound of the parameter for the imposition of Dirichlet type 
boundary conditions. This is related to the fact that the method never uses 
the coecivity of a;i(-, •). 

- Nitsche type boundary terms are added to ah(-,-) in order to ensure con- 
sistency and adjoint consistency, but the penalty is added to the operators 
Sp{-, ■) and Sa{-, ■), allowing for different boundary penalty for the primal and 
the adjoint. As we shall see below for some problems this is the only way to 
make the Nitsche formulation consistent. 

If the primal and the dual problems have a Dirchlet boundary condition on Td this 
is imposed by 

ah{uh,Vh) := a{uh,Vh) - {fiVuh^Vh)^^^ - {^J,yvh,Uh)Y'^ 
and by adding the boundary penalty term 

^DfJ-h^^unVh ds (3.13) 



to Sp{-, ■) and Sa{-, ■) with jjj > 0. In the non- homogeneous case the suitable data is 
added to the right hand side in the standard way. For Neumann conditions on Tn in 
both the primal and the adjoint problems, these are introduced in the standard way 
in a{u, v) with a suitable modification of the right hand side of (|2.ip . No modification 
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is introduced in ah{-, •) but the following penalty is added to Sp(-, ■) and Sq(-, ■) with 
7JV >0 

/ j^hVuii ■ n\7vh • n ds. (3-14) 

"'rjv 

If the CIP method is used the only difference between Sp{-,-) and Sa{-,-) are pos- 
sible difference in the boundary conditions. If the boundary conditions for u are 
non-homogeneous the usual data contributions are introduced in the right hand side 

~Sp{u,Wh). 

As mentioned in the introduction the semi-norm | • |s can be a norm in certain 
situations so that the partial coercivity (|2.10|) implies the well-posedness of the linear 
system ()2.5|) . In the following Proposition we discuss some basic sufficient conditions 
for the matrix to be invertible in the case of piecewise affine approximation spaces. 
For particular case other geometric arguments may prove fruitful as we shall see in 
the second example below. 

Proposition 3.1. The kernel of the linear system defined by (|2.5p with the 
stabilisation (|3.6p has dimension at most 2{d+ 1) for k — I. The system (|2.5|) admits 
a unique solution if the boundary conditions satisfy one of the following conditions: 

1. Two non-parallel polygon sides subject to Dirichlet boundary conditions. 

2. Two non- orthogonal polygon sides subject one to a Dirichlet boundary condi- 
tion and the other a Neumann condition imposed using (j3.14p . 

3. d non-parallel polygon sides subject to Neumann conditions imposed using 
(|3.14p and either I ^ Vh or there exists Vh,Wh G Vu such that ah{l,Vh) ^ 

and Qhiwh,!) ¥"0- 

Proof. It is immediate from (|2.10p that the kernel of the system matrix of (|2.5p 
can not be larger than the sum of the dimensions of the kernels of Sp(-, •) and Sq(-, •)■ 
For Scip{-, ■) and k — 1 the kernel is identified as [Pi(ri)]^, with dimension 2{d -\- 1). 

To prove wellposedness of the linear system it is enough to prove uniqueness, we 
assume that / = and prove that then u/i = z/i = 0. 

If Dirichlet boundary condition is imposed on a boundary then the gradient must 
be zero in the tangential direction to this boundary, since the tangents of two bound- 
aries spans K'' we conclude that / = in (|2.5p implies u^ = due to (|2.10p and 
similarly z^ — and the matrix is invertible. 

In the second case, the function is zero on the Dirichlet boundary and the gradient 
is zero in the tangential directions of the Dirichlet boundary condition eliminating d 
elements in the kernel. The penalty on the Neumann boundary, being non-orthogonal 
to the Dirichlet boundary, cancels the remaining free gradient. The same argument 
leads to both Uh = and Zh = 0. 

For the third case we observe that the term (|3.14p acting on d non-parallel polygon 
sides implies that Vuh = and wellposedness is then immediate by the remaining 
conditions. D 

Remark 3. Observe that the Proposition \377\ holds for any bilinear form a(-,-), 
even strongly degenerate ones. 

4. Applications. We will now give three examples of problems that enter the 
abstract framework. The first two problems we introduce below have well-posed 
primal and adjoint problems so that the above theory applies. For each method we 
will propose a formulation and prove that the relations ()2.7p . (|2.8p hold. We only 
consider the CIP-method in the examples below, in the first example we detail the 
dependence of physical parameters in all norms and coefficients. In the later examples 
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we assume that all physical parameters are unity and do not track the dependence. 
As suggested above we take Try = tt\y = tt^. In each case we will detail the form 
of e{h). In the last case, the elliptic Cauchy problem, the stability properties of the 
problem is strongly depending on the geometry of the problem and the assumption 
of well-posedness does not hold in general. We will nevertheless propose a method 
that satisfies the assumptions of the general theory and then study its performance 
numerically. 

4.1. Nonsymmetric indefinite elliptic problems. Our first examples con- 
sist of a convection-diffusion-reaction problem with non-solenoidal velocity field as 
is the case for reactive transport in compressible flow. We first consider the case of 
homogeneous Dirichlet conditions where the analysis of |23j applies. Then we con- 
sider the case where failure of the coercivity is due also to the boundary condition, 
here we study a convection-diffusion equation with homogeneous Neumann boundary 
conditions. We will detail only how the analysis of this case differs from the Dirichlet 
case. For a detailed analysis of the well-posedness of the continuous problems we refer 
to [m [16] and for a finite element analysis in the case of homogeneous Neumann con- 
ditions to [To]. Recent work on numerical methods for these problems have focused 
on finite volume methods, see [151 111] or hybrid finite element finite volume methods 

m- 

4.1.1. Reactive transport in compressible flow: Dirichlet conditions. In 

combustion problems for example it is important to accurately compute the transport 
of the reacting species in the compressible flow. We suggest a scalar model problem 
of convection-diffusion type with a linear reaction term cu, where the reaction can 
have arbitrary sign. 

(4.1) 



Cu = f in il 




u ~ on dfl. 




.z~P-Vz + cz = g 


in ri 


z = 


on d^ 



The dual adjoint takes the form 

L*z:= -^J^/^z~ 13 ■Vz + cz = g in f^ (4.2) 

z = on do.. 

The variational formulation is obtained by taking V :— H}^{Vt) and introducing the 
bilinear form 

a{u, v) :— (/iVu, Vv) + (V ■ (/3u) -I- cu, v). 

We assume that f,g & L'^{^), that both (|4.1|) and (|4.2|) are well posed in Hq{^) by 
the Fredholm alternative and that the smoothing property p.3|) holds. See [14] for an 
analysis of existence and uniqueness under weaker regularity assumptions on /3 and 
c, with c > 0. The below analysis can also be carried out assuming less regularity, 
but the constraints on the mesh-size for the error estimate to hold will be stronger. 
Recall that the constants in the estimate ()2.3p also depends on the regularity of the 
coefficients. 

The discrete form of the bilinear form is given by 

a/i(w/i, Vh) := a{uh,Vh) - (Vu/i ■ n, Vh)Qf^ ~ {Vvh ■ n, uh) q,^ - {{P ■ n)_Uh,Vh)g^ (4.3) 

where {(3 ■ n)± := ^(J3 ■ n±\/3 ■ n\). The stabilisation is chosen as 

Spiuh, Vh) ■■= Scipiuh, Vh) + S'bci'^h, Vh) (4.4) 
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and 

Sa{zh,Vh) := Scip{zh,Vh) + S^^{zh,Vh) (4.5) 

where 71. f '^ (M+||/3h''T-F||oo,_F^F) and 72, f '^ /i in p.6p with /3fi the nodal interpolant 
of /3 and 

s^c(^''''"'i) := {f^h^^Xh,Vh)g^ + {\{P ■n)±\xh,Vh)gQ. (4.6) 

Proposition 4.1. (Existence of discrete solutions) Let k ~ 1. Then the for- 
mulation (j2.5p with the bilinear form (I4.3p and f/ie stabilisation (j4.4p - (|4.5p admits a 
unique solution u^ 6 V/j. 

Proof. Immediate by Proposition 13. II D 
It is well known that the bilinear form (|4.3p satisfies the consistency relations (|2.7p and 
(|2.8p and that the stabilisation (|4.4p - (|4.5p satisfies the consistency properties (|2.13p . 
()2.14p and p.lSp . Now we define the norms by 

ll«ll+ := \M* := ll/iU-^wll^r^^^ + IKCpe + Cosc + Co)«|| + ll/i^A^^ Vi;||an + |k^f Han- 

HereCpe := (cr^ft-^^ +Ai^/i"^), with cr := ||/3|| 1,00(0), Cosc := /i"^(|^|w2,~ + |c|vi/i.-)^ 
and Co = (||V/3||l°= + ||c||l°°)^- Also define 

II^IU := U^h/^vWh + ||(7~U^/3 • Vv|| + |||c|5i,|| 

+ Wii^h^lVv ■ nFl||.F„,, + W~^h--2 + a5)t;||ao + \v\s, + e(/i)||w|!. 
It is straightforward to show that 



\\u ~ ■nvu\\+ + \\u-t:vu\\c < {(pe + (osc + Co)h^^^\u\Hk+i 



(n) 



and similarly for ||m — tt^vmH*. 

It then only remains to prove the continuities (|2.1ip and (|2.12p to conclude that 
the Theorem!^ holds. 

Proposition 4.2. The bilinear form (|4.3p satisfies the continuities (|2.1ip and 
(Pl^ with 

e{h) -/i5(|/3|,^2,oo + |c|vi/i,-)^ 



Proof. First we consider p. lip . After an integration by parts in a(-, •) we have 

ah{u-'Kvu,Xh)^ ^ {u-'Kvu,l^JSIxh■npl)p + ^{u-■nvu,C*Xh)K 

- (u - TTyu, (^ • n)+a;/i)gs-j = I + II + III. 

Considering / — /// we find using Cauchy-Schwarz inequality 

I + III < \\u- nvu\\ + \xh\s^- 

For //, using the discrete interpolation results p.lOp . the discrete commutator prop- 
erty, see |3] , and standard approximation followed by an inverse inequality in the last 
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term 

K 
+ (m - TTyM, CXji - ihicXh)) + (U - TTyti, (;3 - 4;3) • Vx/i) 

< Cos\\u - TTvu\\ + {\xh\s^ +Ci/l^(|/3|H'2,oo + |c|H'l,oo)3||a;,,||). 

The second continuity follows in a similar fashion, 

ah{u-Uh,y~Trwy)^{C{u-Uh),y~TTwy)h+ ^ ([a*Vu/i • n^l, J/ - tth/?/)^ 

+ ((/3 • 'n)^Uh, y - T^wy)sn + if^^iv - '^wv) ■ n, Uh)g^ 
= 1 + 11 + 111 + IV. 

Considering first the term / we get, with ^h = Tryu — Uh 

I = {C{u - TTyu) + CS,h, y - TTwy) 
^ \\y - T^wyhiWu - TTyuWc + |6i|sp +h'^{\P\w^.--{n) + \c\w^.'^{n})^\\u - Uh\\), 
where we used once again the inequalities 

ifiA^h + Ph^ih,y-T^wy) < Cos\ih\s^\\y - T^wyW*, 

{{f3 -ihP) ■V^h,y-T^wy) < h\/3\w2.oa{\\u-TTvu\\ + \\u - Uh\\)\\y - irwyW 

and 

((V • /3 + c)^h,y - TTwy) = ((V • /3 + c)^h ~ Jh((V • /? + c)^), y ~ ^wy) 

< h{\(3\w2.oa + |c|wi.oo)(||w - TTyuW + \\u - uh\\)\\y - TTwyW- 

In the last inequality we once again applied then discrete commutator property. For 
the second term and third terms we have using the Cauchy-Schwarz inequality, re- 
calling the form of the boundary penalty term, 

ii + iii + iv < \uh\sjy - T^wyW*- 

Hence, by adding and subtracting nyu in the \uh\sp contribution and using approx- 
imation we conclude that the claim holds with e{h) ~ /12 (|/3|^y2,oo -|- |c|vi/i.°°)^- D 



Remark 4- Note that if the physical parameters are constant, then the analy- 
sis holds without restrictions on the mesh size in contrast to the standard Galerkin 
analysis of fl^ . 

Remark 5. There is some flexibility in the choice of norms in the analysis above. 
Here we have tried to mimic the analysis of a stabilised method that is robust for 
high Peclet number, however this estimate is not robust, even if we assume that u 
is as smooth as we like. This is due to the use of a Nitsche type trick where the 
actual regularity constant of the estimate (|2.3p appears. If the dependence on physical 
parameters in Theorem \2.2\ is taken into account the upper bound typically reads for 
k = l, 

IW - Uh\\ < Ca.niCPe + Cosc + Co)^^" I"|ff2 (n) • 
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Assuming that all physical parameters except /i are 0{1) we get 

\\u -Uh\\ ^ ca,n ( r + ri ) h*\u\H2(n) 

Here the constant Ca,n typically is proportional to some negative power of fj,, making 
the estimate valid only for moderate Peclet numbers. If we assume that Ca.n = 0{^^^) 
we see that the quasi optimal convergence of order h^ is obtained when h"^ < fi. A 
more precise estimate for the hyperbolic regime is the subject of the second part of this 
work 1^. 

4.1.2. Transport in compressible flow: pure Neumann conditions. We 

will now consider the convection-diffusion equation with homogeneous Neumann con- 
ditions. The main difficulty in this problem compared to the previous one is that due 
to the homogeneous Neumann condition, the primal and dual problems have different 
boundary conditions. The non-solenoidal /3 imposes special compatibility conditions 
on g leading to complications in the finite element analysis and additional stability 
issues for the discrete solution. For this example we will assume that all physical 
parameters are order one. After having presented the problem and the method we 
propose, we first show that the the discrete problem is well-posed for all mesh-sizes 
when piecewise affine approximation is used. Then we prove that the assumpitions 
of Lemma 12.11 and Theorem 12.21 are satisfied. Optimal error estimates for the prob- 
lem similar to that above are obtained after accounting for some minor modifications 
needed to accomodate the compatibility conditions particular to this problem. The 
problem reads 

-Au -t- V • {I3u) = f inn (4.7) 

— Vu ■ n + /3 ■ nu = on dfl. 

The dual adjoint problem is formally written 

-Az - p-Vz^g infl (4.8) 

— Vz ■ n = on dil. 

We assume that the following compatibility conditions hold 

/ / da: = 0, I gmdx^Q (4.9) 

where m G H'^iyi), m > is the unique solution to the homogeneous form of the 
primal problem 

-Am 4- V • (/?m) == in fi (4.10) 

— Vm ■ n -\- P ■ nm — on dfl. 



under the additional constraint 



\n\-^ mdx 
Jfi 



Then the problems (|4.7|) and ()4.8p are both well-posed by the Fredholm alternative. 
Since we assume that the regularity estimate (12. 3p holds, m G C°{Cl) and sup^g^ ^'^ =• 
M G R+ and since m > we may introduce JTimm '■= infn m > (see |10jV 
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The problem is cast in the form (|2.ip by setting V :— iJ^(il) n LQ{fl), where 
io(^) denotes the set of functions with global average zero, and by taking 

a{u, v) := (Vm, Vv) - [u, i3 ■ Vv). 

The finite element method (|2.5|) is obtained by setting Vh :— Wh '■— X^ n Lq{VL) 

ah{uh,Vh) -^ a{uh,Vh) (4-11) 

and the stabilisation operators 

s^(-, •) := Sip{-, ■) + Sbc,x{'^ •)) with X ^ a,p. (4-12) 

Scipi:, ■) is given by (|3.6p with "/i,p := 1, i = f , 2. The boundary operators finally are 
defined by 

Sbc,p{uh,Vh) ■■= I h{Vuh -n- 13 ■ nuh){yvh ■ n - 15 ■ nvh) ds 
Jo. 

for fc > 2 and 

Sbc,piuh,Vh) := / h{\/uh ■ n - (ihP) ■ nUh){Vvh ■ n - [ihjS) ■ nv^) As 
Jn 

for fc = f , Sbc,a{', ■) finally is given by equation (|3.f4p . with jn = 1. The boundary 
stabilisation operator for fc = 1 is only weakly consistent. It is straightforward to show 
that the inconsistency introduced by replacing (3 by i^P is compatible with (|2.f3l) . We 
omit the details here, but similar arguments are used below to prove the continuity 

Proposition 4.3. (Existence of discrete solution) Assume k = 1 in the definition 
ofVh- Then there exists a unique solution Uh to the discrete problem ()2.5p . 
Proof. As before we assume / = and observe that 

Sp{uh,Uh) + Sa{zh,Zh) = 0. 

This implies Zh,Uh G Pi(fi). Since ||Vzft, • n||asi = and Sav{zh,Zh) — 0, we conclude 
that Zh = 0- For Uh there holds 

\\Vuh ■ n + {ihP ■ n)uh\\dn = 0. 

Since Vuh ■ n is constant on every polyhedral side F of ri so is (j/j/3 • n)uh- But since 
(i/j/? • n)uh\r G IP'2(F) we conclude that both ihl3 and Uh must be constant. Since this 
is true for all faces F of J7, Uh is a constant globally. We conclude by recalling that 
zero average was inposed on the approximation space. D 

In case fc > 2, we let the norms || • || + , || • ||* be defined by ^^ and || • 1|£ by (|XT^ . 
When fc = 1 we let the norm || • ||+ be defined by (|3.9p . but define 



\\v\U:=\\h-\\\ + \\h-^v\\^ + eih)\\v\\ 

and 

\\v\\c ■■= \\Cv\\h + Wh^Vv ■ npE^.^t + l|/i*Vw • nWgn + {e{h) + (T)\\h^v\\9n + e{h)\\v\\. 

For the projection operators Try and t^w we once again choose the L'^-projection. 
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Proposition 4.4. Assume 

e{h)^h^\/3\w2,o.. 

Then the bilinear form (|4.1ip satisfies the continuities (|2.1ip and ()2.12p . 
Proof. As before we integrate by parts in a/j(-, •) to obtain 

ah{u-TTvu,Xh) = ^ {u-Trvu,l\7xh-nFi)p + ^^{u- nvu,C*Xh)K 
FeJ^M K 

+ {u — TTyu, Vx/i • n)gj2 ~ I + II + III- 

The treatment of terms / and // are identical to the Dirichlet case. Term /// is 
bounded using Cauchy-Schwarz inequality, recalling that the Neumann condition is 
penahsed in Sa{-,-) 

III < \\u-Trvu\\ + \xh\s^- 
The second continuity follows in a similar fashion. We write 

ah{u - Uh, y - TTwy) = ^ mu - un), y - T^wy)K 

+ X! {\^^h-nFl,y-i^wy)F 

- (y - TTwy, ^Uh -n- {13 ■ n)uh)gQ 
= 1 + 11 + 111 

and observe that the treatment of terms / and // is analoguous with the Dirichlet 
case. Recalling that Vw^ ■ n — {j3 ■ n)uh is penalised in Sp(-, •), we may conclude as 
before using a Cauchy-Schwarz inequality when k > 2. 

Ill < \\u-Trvu\\^\uh\s^. 

For the case fc = 1 we must take care to handle the lack of consistency. Therefore we 
add and subtract ifif3 and use the boundary condition on u to get 

/// = {y - TTwy, ^{uh -u) -n- [ihfi ■ n){uh - u))q^ 
+ {y~ TTwy, {ihP - (3) ■ n{uh ~ u))g^^ . 

First we add and subtract wvu, so that u — Uh= u — irvu+ih, ^h '■— ttvu — uh and split 
the scalar products with Cauchy-Schwarz inequality. Then we apply a trace inequality 
in the boundary norm on ^^, and use the definition of the norms, in particular that 
ll^~^(2/ ~ 7''TV2/)||ao < \\y — TTwyW* and apply approximation in _L°° for j3 to obtain 

III ^ \\y - T^wy\\*{\\u - t:vu\\c + l^/ilsp) + \\y - TTwy\\d<:ih'^\P\w^''^\\ih\\- 
Using a triangular inequality in \\£,h\\, after having added and subtracted u, we obtain 

III ^ \\y - T^wyhiWu - irvu\\c + |6i|sp) + \\y - ■TTwy\\dnh^\P\w^MW ~ '^hW- 
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Using the definition of the norms we see that the last term in the right hand side may 
be written 

\\y - T^wyhnh^lPlw.-^Wu - Uh\\ < \\y - Trwy\\*h'^\l3\w^.-^\\u - Uh\\ 

~ \\y - '!Twy\\*e{h)\\u - Uh\\. 

The proof is complete. D 

Remark 6. It is straightforward to verify that Lemma \2. 1\ holds. The assumptions 
of Theorem \2.2\. however still are not satisfied since we want to use the solutions of the 
problems C*ip = u^Uh and Cip = Zh, but the solution Lp will in general not exist since 
u — Uh does not satisfy the second compatibility condition of equation (|4.9p . Instead 
we will use m, the solution of (J4.10I) as weight, as suggested in ^1 OJ . and solve the 
wellposed problem 

C*ip = (u — uii)/m. 

We may then write 

|(u - -u/i)m"2 1|2 _|_ ||zft|p = (-U - Uh,C*ip) + {C'4>,Zh) 

and proceed as in Theorem \2.2l now using the stability estimate 

IflHHn) < Ca^Wiu ~ Uh)m-^\ < Ca,nM™„||(w - Uh)m-^/^\\ 

to obtain 

\\{u~Uh)m-^\ + \\zh\\ < h''+^\u\H^^+l(n)■ 

For an estimate in the unweighted L^-norm we observe that 

M-^||u-u,,|| < \\{u~Uh)m-i\\. 

Convergence follows by Lemma \2.1\ and the modified Theorem \2.2\ Observe that the 
constants in e(h) now depends on the (unknown) minimum value of m. 

Remark 7. In practice the zero average condition can be imposed using Lagrange 
multipliers. The above analysis holds for that case after minor modifications. 

4.2. The Cauchy problem. We consider the case of a Helmholtz-type problem 
where both the solution itself and its normal gradient is specified on one portion of 
the domain and the other portion is free. We let Ty and T\y be connected subsets 
of dVl such that dVl :— Ty U Tw and Ty n Tw = 0. We will consider the problem, 

-Aw + Ku = /inrj (4.13) 

Vm • n = m = on Ty, 



with dual problem 



-Az + Kz = g in fl (4-14) 

Vz • n = z = on Tw- 



The weak formulations (|2.1I) and (12.21) are obtained by setting 

V ■.= {ve H\n) : v\r^r = 0} 



19 



and 

W := {v e H\n) : v\r„ = 0} 
and defining 

a{u, v) := (V-u, Vv) + k{u, v), \^u £ V, v £ W. 

Note that both symmetry and Gardings inequality fail in this case because the func- 
tions in the bilinear form have to be taken in different spaces and hence the choice 
V = u is prohibited. 

To design a suitable discrete formulation (|2.5p for this problem we generalise the 
ideas of the Nitsche type weak imposition of boundary condition. Observe that in 
this case boundary conditions imposed using penalty in the standard fashion can not 
be consistent for both the primal and the adjoint problems, since the primal and dual 
solution are zero on different parts of the boundary. It is therefore important in this 
case that two stabilisation operators are used, one for the primal and one for the 
adjoint. We propose the bilinear form ah{-, ■) ■ Vh x Vh -> R defined as 

ah{uh,Vh) ■^a{uh,Vh) - {Vv^ ■n,u,i)^^ - (Vu,, ■n,Vh)^^^ (4.15) 

and for the stabilisation we use 

Sx{uh,Vh) := Scipiuh,Vh) + Sbc,xiuh,Vh), x = a,p (4.16) 

where Scip{-, ■) is given by (|3.6p . with ^p_i :— 1^ i — 1, 2, and 

Sbc.x{uh,Vh) -^ / {h^^UhVh + hVuh-nVvfi-n) ds, 
J X 

where X = Ty for x = p and X = Tw for x — a. If some part of the boundary is 
equipped with Dirichlet or Neumann boundary conditions this is imposed as described 
in Section [331 

Proposition 4.5. (Existence of discrete solution for k — 1) Define (|2.5|) by the 
bilinear forms (|4.15p and ()4.16|) . Let k — \ inVh- Then there exists a unique solution 

iuh,zh) e[Vh]^ to dm. 

Proof. Let / = 0, by (|2.10p there holds, u/i, z^ G Pi(r2) and Uh\rv — ^^h ■ n\rv = 
as well as ^^/iItm,' — Vz/i-nlri,,- = 0, by which we conclude that the matrix is invertible 
using case (2) of Proposition [3Jl D 

For the error analysis we once again choose the interpolants Try and ttw to be the 
standard L^-projection ttl. We will now prove that the assumptions ()2.7p - p.8|) and 
(P1T|) - (PT^ are satisfied. 

Lemma 4.6. (Consistency of bilinear form) The bilinear form (J4.15P satisfies 
((2Jl) and (gj]). 

Proof. By an integration by parts we see that for u solution of ()4.13p 

(-A?i + Ku,v + Vh) = (Vu, V(w + Vh)) - (Vu ■n,v + Vh)j. - (V(w + Vh) ■ n, u)p 



since Vm n=o on Fv since u=o on Tv 

+ {ku,v + Vh) = ah{u,v + Vh). 

Similarly for z solution of (|4.14l) consistency follows by observing that 

(u + w/i,-Az) = (V(w + w/i), Vz) - (Vz •n,w + t;/i)p - (V(u + «;>)• n, z)p 



since V2n=o on r^ since z=o on Tw 
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D 

We define the norms || • I|+ and I| • ||, by 

\\v\\+ := Wh-hy^^^ + \\h-'v\\ + \\h-h\\r,y + \\h-^Wv ■ n\\r^, 

\\v\U ■■= \\h-h\\^^^, + \\h-^v\\ + \\h-K\\r,. + Wh^Vv ■ n\\r^. 
and 

\\v\\c := \\hCvU + WhH^v ■ M|b.„, + \\h-h\\r, + \\h-^Vv ■ n\\r,. 

It is straightforward to show (|2.13p and (I2.14p . 

Proposition 4.7. For a/i(-,-) defined by (|4.15p . the continuities (|2.1ip and 
(Pl^ hold with e{h) = 0. 

Proof. We proceed as before using an integration by parts in (I4.15P to obtain 

ah{v - TlvV, Xh) = X! ^^ ^ '^^^' l^^h ■ nFi)p +^{v - Try W, -Axh + KXh)K 

+ {v ~ TTyv, Vxh ■ n)p . — (V(w — TTyv) ■ n, xii)-p , = I + II + III + IV. 

The first sum I is upper bounded as before using the Cauchy-Schwarz inequahty and 
for the second sum, we use the orthogonahty of the L^-projection, {v — Try w, KXh) = 
and the discrete interpolation inequality p.lOp leading to 

I + II ^ \\u-TTvu\\ + \xh\sa- 

For the terms /// and IV we note that by the definition of || • ||+ and | • l^^ there also 
holds 

III + IV < ||u-7ryu|| + |a;,,|s„. 

This ends the proof of (|2.1ip . The proof of p.l2p is similar. Using integration by 
parts in the other direction we have 

Ohiu ~uh,y - TTwy) = ^{-M^ - "'i) + «(w ~ w/i), y ~ Trwy)K 

K 

+ X! (iV-u/i ■ni;^],?;-7ri4/?;)^ + (u-'U/i,V(y-7rwy) ■n)p^ 

+ (V(u ~Uh)-n,y^ t^wv)t^- =I + II + III + IV. 

Using the same arguments as before, adding and subtracting TTyu in all the terms in 
the left slot we have for the term /, using i^^ = TTyu — u^, 

I = (-A(u - TTvu) + k{u - TTyu), y - TTwy)h - (A^/i - IosA^h,y - Trwy)h 
< {\\u - TTvuWc + \£.h\sp)\\y - T^wyW*- 

By the definition of the stabilisation operator and the fact that u = Vu ■ n = on Fy 
we may once again add and subtract TTyu in the terms // and /// to obtain 

// + /// + IV= ilVuh ■nFJ,y- T^wy)^^^^ 

- {uh, V(y - TTwy) ■ n)p^^ - {Vuh ■n,y- Trwy)r^ 

< \\y - T^wyhilW - TTvuWc + IC/ilsJ- 
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By which we conclude. D 

Corollary 4.8. Assume that the problems (j4.13p and (j4.14p adm,it unique so- 
lutions for which (|2.3p holds. Then the conclusions of Theroem \2.'2i holds for (|2.5p 
defined by Vh, (|4A5| and (|4T6)) . 

Proof. In Lemma l476l we verified the consistencies (|2.7I) and (|2.8p . In Proposition 
14.71 we verified the continuities (|2.1ip and ()2.12p . It is straightforward to verify that 
(|2?T3ll - (f214l) hold for ny = ttw = t^l and Sp{-, •), Sa(-, •) defined by (j4?T6l) under the 
assumptions on the mesh and the regularity assumptions on the solution. D 

Remark 8. Admittedly Corollary |4.<^| is of purely academic interest since the 
Cauchy problem under consideration in general is ill-posed, with very weak stability 
properties. As we shall see in the numerical section the method nevertheless returns 
useful approximations. An example of a sufficient condition for Theorem \2.2\ to result 
in a convergence estimate is that there exists M S M"*" and s £ R, with s > — fc/2 such 
that ||Vz|| < Mh^ , for all u — u^,. The expected convergence rate in that case would 
be 

\\u-Uh\\<h''/^+'\u\H.+^n), (4.17) 

however since ||Vz|| can be expected to vary significantly for different realisations Uh 
the convergence must be expected to be very irregular, even if (|4.17p holds. Unfortu- 
nately, no such stability estimates are known for the Cauchy problem, we refer to ^ 
for conditional stability estimates for the problem (|4.13p in general Lipschitz domains 
and to \n\ \2i^ and I21'l for other work on finite element methods on the Cauchy 
problem and some stability results under special geometrical assumptions. 

5. Numerical investigations. We will present numerical examples of conver- 
gence for a smooth exact solution of the applications given above. For the computa- 
tions we have used FreeFEM-|— 1-, [TH]. All problems will be set in n := (0, 1) x (0, 1). 
We consider an example given in [11] where, in (|4.ip . the physical parameters are 
chosen as^ = 1, c = 0, 

/3:=-100f ^ + ^ 
\ y-x 

and the exact solution is given by 

uix,y)^30x{l-x)yil-y). (5.1) 

This function satisfies homogeneous Dirichlet boundary conditions and has ||m|| = 
1. Note that ||/3||ioo = 200 and V • /3 = —200, making the problem strongly non- 
coercive with a medium high Peclet number. The right hand side is then chosen 
as Cu and in the case of (non-homogeneous) Neumann conditions, a suitable right 
hand side is introduced to make the boundary penalty term consistent. We use 
unstructured meshes with 2-^^ elements on each side, N = 3, . . . , 8 and, drawing 
on our previous experience of the CIP-method we fix the stabilisation parameters 
to be 7i.F = 0.01 for piecewise affine approximation and ^i^p = 0.001, i = 1,2 for 
piecewise quadratic approximation. The boundary penalty parameter is chosen to be 
7f,c ~ 10 for both cases and both for Dirichlet and Neumann penalty terms. Let us 
remark that in particular for the ill-posed Cauchy problem, an optimal choice of the 
stabilisation parameter can make a big impact on the error on a fixed mesh, but does 
not appear to influence the convergence behavior. For each example we plot the error 
quantities estimated in Lemma [2.1l and Theorem l2.2l When apropriate we indicate the 
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Fig. 5.1. Left: example of unstructured mesh, N = 5. Right: plot of the velocity vector field. 



N 


1 "" ~ ""/i 1 


ll^/ill 


\uh\Sp + \zh\Sc, 


3 


0.038 (-) 


0.024 


0.57 


4 


0.012 (1.7) 


0.0017 


0.24 


5 


0.0024 (2.3) 


0.00043 


0.11 


6 


0.00043 (2.5) 


0.00012 


0.052 


7 


0.00010 (2.1) 


2.5E-05 


0.025 


8 


2.3E-05 (2.1) 


5.3E-06 


0.012 



Table 5.1 

Convergence orders of estimated quantities for the Dirichlet problem approximated using piece- 
wise affine elements 



experimental convergence order in parenthesis. We report the computational mesh 
for iV = 5 and the velocity field (3 in figure [5l The optimal convergence rate for the 
stabilizing terms given in Lemma l2. H is verified in all the numerical examples. 



5.1. Dirichlet boundary conditions. In table ISTTJ we show the result of the 
computation when Dirichlet boundary conditions are applied and piecewise afHne 
approximation is used on a sequence of unstructured meshes. We observe that the 

3 

solution exhibits the preasymptotic convergence rate h^ under one refinement before 
achieving the full second order convergence rate in L^. 

In table 15.21 similar data for second order polynomials are presented. Here the 
asymptotic regime with full convergence is obtained from the first refinement. 

5.2. Neumann boundary conditions. We consider the same differential op- 
erator, but with (non-homogeneous) Neumann-boundary conditions. This is exactly 
the problem considered in |11) . The average values of the approximate solutions have 
been imposed using Lagrange multipliers. In tables 15.31 and 15.41 we observe optimal 
convergence rates once again as predicted by theory. Observe that in the case of 
piecewise affine approximation the dual solution Zh comes into the asymptotic regime 
only on the finer meshes. 

5.3. A Cauchy problem. Since we have no complete theory for the ill-posed 
Cauchy problem we will proceed with a more thorough numerical investigation. First 
we consider the Cauchy problem obtained by taking k = in (I4.13p . Then we consider 
a Cauchy problem using the convection-diffusion operator of (j4.7p in two different 
boundary configurations. For all test cases we use the exact solution (|5.ip and the 
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N 


\\u — UhW 


\\zh\\ 


\uh\sp + \zh\s^ 


3 


0.0014 (-) 


0.00041 


0.024 


4 


0.00012 (3.5) 


4.6E-05 


0.0044 


5 


8.8E-06 (3.8) 


4.6Ee-06 


0.00081 


6 


8.0E-07 (3.5) 


6.6E-07 


0.00017 


7 


8.3E-08 (3.3) 


8.2E-08 


3.7E-05 



Table 5.2 

Convergence orders of estimated quantities for the Dirichlet problem approximated using piece- 
wise quadratic elements 



N 


\\u — Uh\\ 


W^hW 


\Uh\Sp + \Zh\Sa 


3 


0.028 (-) 


0.028 (-) 


0.82 


4 


0.0066 (2.1) 


0.016 (0.8) 


0.32 


5 


0.0016 (2.0) 


0.0058 (1.5) 


0.13 


6 


0.00039 (2.0) 


0.0015 (2.0) 


0.060 


7 


9.7E-05 (2.0) 


0.00031 (2.3) 


0.028 


8 


2.3E-05 (2.1) 


6.5E-05 (2.3) 


0.013 



Convergence orders of estimated 
wise affine elements 



Table 5.3 

for the Neumann problem approximated using piece- 



stabilisation parameters given above. We present the data for the quantities estimated 
in Lemma 12.11 and Theorem 12.21 but also the error in the total diffusive flux in the 
discrete H^^'^{dfl) norm on the boundary. 



\\V{u-Uh) -nll^i ^ 



on 



ft-(V(u — Ufi) ■ n) ds. 



on 



5.3.1. Poisson's equation. Here we consider the problem with k = in (|4.13p . 
We impose the Cauchy data, i.e. both Dirichlet and Neumann data, on boundaries 
a; = 0, 0<y<l and y = 1, < x < 1. In table [5?2] we show the obtained errors 
when piecewise affine approximation is used and in table 15.21 the results for piecewise 
quadratic approximation. 

First note that in both cases one observes the optimal convergence of the sta- 
bilisation terms predicted by Lemma 12.11 For the L^-norm of the error we observe 
experimental convergence orders /i" with typically a ^ 0.25 for piecewise affine ap- 
proximation and a ^ 0.5 for quadratic approxmation. Higher convergence orders were 
obtained in both cases for the normal diffusive flux. In figure l^?^ we present a study of 
the L^-norm error under variation of the stabilisation parameter. The computations 
are made on one mesh, with 32 elements per side and the Cauchy problem is solved 
with k = 1,2 and different values for 'jf,i = 7_f,2 with 76c = 10 fixed. The level of 10% 
relative error is indicated by the horizontal dotted line. Observe that the robustness 
with respect to stabilisation parameters is much better for quadratic approximation. 
Indeed the 10% error level is met for all parameter values ji^p G [2.0E — 5, 1], whereas 
in the case of piecewise affine approximation one has to take 7i_f G [0.003,0.05] ap- 
proximately. Similar results for the boundary penalty parameter not reported here 
showed that the method was even more robust under perturbations of jbc- In the 
right plot of figure [?31 we present the contour plot of the error u — Uh and the contour 
plot of Zh- In both cases the error is concentrated on the boundary without boundary 
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N 


llw-w/ill 


\\Zh\\ 


\uh\sj, + \zh\sc. 


3 


0.00061 (-) 


0.0020 (-) 


0.030 


4 


6.6E-05 (3.2) 


0.00040 (2.3) 


0.0054 


5 


6.5E-06 (3.3) 


2.5E-05 (4.0) 


0.00099 


6 


7.1E-07 (3.2) 


1.7E-06 (3.9) 


0.00020 


7 


7.9E-08 (3.2) 


1.4E-07 (3.6) 


4.2E-05 



Table 5.4 

Convergence orders of estimated quantities for the Neumann problem approximated using piece- 
wise quadratic elements 



N 


\\U-Uh\\ 


lk/.|| 


\Uh\s^ + \Zh\Sa 


\\V{u-Uh)-n\\_i,^Q^ 


3 


0.070 (-) 


0.59 (-) 


2.0 (-) 


2.7 (-) 


4 


0.074 (-) 


0.42 (0.49) 


0.79 (1.3) 


1.3 (1.1) 


5 


0.037 (1.0) 


0.30 (0.49) 


0.30 (1.4) 


0.75 (0.80) 


6 


0.029 (0.35) 


0.26 (0.2) 


0.13 (1.2) 


0.51 (0.56) 


7 


0.024 (0.27) 


0.20 (0.37) 


0.054 (1.3) 


0.33 (0.62) 


8 


0.020 (0.26) 


0.16 (0.32) 


0.022 (1.3) 


0.21 (0.65) 



Table 5.5 

Convergence orders of estimated quantities for the Poisson Cauchy problem approximated using 
pieceuiise affine elements 



condition. 

5.3.2. The non coercive convection diffusion equation. As a last example 
we consider the Cauchy problem using the non-coercive convection-diffusion operator 
()3.ip . The stability of the problem depends strongly on where the boundary conditions 
are imposed in relation to the inflow and outflow boundaries. To illustrate this we 
propose two configurations. Recalling the right plot of figure [5] we observe that the 
flow enters along the boundaries y = 0, y = 1 and x = 1 and exits on the boundary 
x = 0. Note that the strongest inflow takes place on y = and a; = 1, the flow being 
close to parallel to the boundary in the right half of the segment y = 1. We propose 
the two different Cauchy problem configurations: 
Case 1. We impose Dirichlet and Neumann data on the two mixed boundaries a; = 

and y = 1. 
Case 2. We impose Dirichlet and Neumann data on the two inflow boundaries y = 

and X = \. 
In the flrst case the outflow portion or the inflow portion of every streamline are 
include in the Cauchy boundary whereas in the second case the main part of the inflow 
boundary is included. This highlights two different difficulties for Cauchy problems 
for the convection-diffusion operator, in case 1 we must solve the problem backward 
along the characteristics, essentially solving a backward heat equation, whereas in 
case two the crosswind diffusion must reconstruct missing boundary data. 

In tables 15.71 - 15.101 we report the results on the same sequence of unstructured 
meshes used in the previous examples for piecewise affine and piecewise quadratic 
approximations and the two problem configurations. First note that in all cases the 
result of Lemma EH] holds as expected. Otherwise the method behaves very differently 
in the two cases. For case 1 we observe better convergence orders than in the case 
of the pure Poisson problem, typically h'^ for affine elements and h for quadratic 
elements in the L^-norm. Even higher orders are obtained for the global diffusive 
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N 


\\u-Uh\\ 


W^hW 


\uh\s^ + \zh\s^ 


\\V{u~Uh) ■n\\_if^Q^^ 


3 


0.031 (-) 


0.062 (-) 


0.073 (-) 


0.92 (-) 


4 


0.022 (0.49) 


0.025 (1.3) 


0.014 (2.4) 


0.48 (0.94) 


5 


0.013 (0.76) 


0.014 (0.84) 


0.0025 (2.5) 


0.24 (1.0) 


6 


0.0088 (0.56) 


0.011 (0.35) 


0.00047 (2.4) 


0.13 (0.88) 


7 


0.0069 (0.35) 


0.0067 (0.72) 


8.8E-05 (2.8) 


0.080 (0.70) 



Table 5.6 

Convergence orders of estimated quantities for the Poisson Cauchy problem approximated using 
piecewise quadratic elements, ■y = 0.001, ■yi,^ = 10 




0.00001 0.0001 0.001 0.01 

penalty parameter 



Fig. 5.2. Study of the L^ -norm error under variation of the stabilisation parameter, circles: 
affine elements, squares: quadratic elements 



flux in the discrete H^2 norm. The dual variable Zh on the other hand has very 
poor convergence, although it is quite small on all meshes in the case of quadratic 
approximation. Case 2 (control on main part of the inflow) is clearly much more 
difficult. Convergence orders for both the affine case and the quadratic case are poor 
(around ~ /is) and uneven. The diffusive flux still converges approximately as h^ in 
both cases. We conclude that the Cauchy convection-diffusion problem is much less 
ill-posed if for each streamline either the inflow part or the outflow part lies in the 
controlled zone. The fact that we in case 2 controls more of the inflow boundary is 
unimportant compared to the fact that both the inflow and the outflow are unknown 
in the boundary portion around the corner (0, 1). 

6. Concluding remarks. We have proposed a framework for the design of sta- 
bilised flnitc clement methods for non-coercive and non-symmetric problems. The 
fundamental idea is to use an optimisation framework to select the discrete solution 
on each mesh. This also opens new venues for inverse problems or boundary con- 
trol problems, where Tichonov regularisation can been introduced in the form of a 
stabilisation operator with optimal weak consistency properties, eliminating the need 
to match penalty parameter and mesh size to obtain optimal order. The method 
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Fig. 5.3. 
(right plot). 



Contour plots of the error u — ui^ (left plot) and the error in the dual variable zj^ 



N 


\\u~Uh\\ 


\w\\ 


\Uh\Sp + \Zh\Sa. 


\\V{u-Uh) ■n\\_^j^Q^-^ 


3 


0.032 H 


0.044 (-) 


1.6 (-) 


0.35 (-) 


4 


0.010 (1.7) 


0.020 (1.1) 


0.61 (1.4) 


0.13 (1.4) 


5 


0.0045 (1.2) 


0.034 (-) 


0.24 (1.3) 


0.048 (1.4) 


6 


0.0035 (0.36) 


0.052 (-) 


0.10 (1.3) 


0.018 (1.4) 


7 


0.0039 (-) 


0.056 (-) 


0.045 (1.2) 


0.0074 (1.3) 


8 


0.0026 (0.58) 


0.059 (-) 


0.020 (1.2) 


0.0031 (1.3) 



Table 5.7 
Convergence orders of estimated quantities for the convection-diffusion Cauchy problem ap- 
proximated using piecewise affine elements (case 1) 



N 


\\U-Uh\\ 


W^hW 


\uh\sp + \zh\s^ 


\\\/{u-Uh) ■n\\_if^g^ 


3 


0.13 (-) 


0.032 (-) 


1.74 (-) 


0.44 (-) 


4 


0.097 (0.42) 


0.012 (1.4) 


0.63 (1.5) 


0.23 (0.94) 


5 


0.075 (0.37) 


0.010 (0.26) 


0.24 (1.4) 


0.11 (1.1) 


6 


0.067 (0.16) 


0.010 (-) 


0.10 (1.3) 


0.070 (0.65) 


7 


0.063 (0.089) 


0.0097 (0.044) 


0.043 (1.2) 


0.047 (0.57) 


8 


0.056 (0.17) 


0.0082 (0.24) 


0.018 (1.3) 


0.030 (0.65) 



Table 5.8 

Convergence orders of estimated quantities for the convection- 
proximated using piecewise affine elements (case 2) 



fusion Cauchy problem ap- 



N 


\\U-Uh\\ 


W^hW 


\Uh\Sp + \Zh\Sa 


\\'^{u-uh)-n\\_i^f„dn 


3 


0.0022 ( ) 


0.0037 (-) 


0.096 (-) 


0.033 (-) 


4 


0.00054 (2.0) 


0.00089 (2.1) 


0.020 (2.3) 


0.0091 (1.9) 


5 


0.00024(1.2) 


0.0013 (-) 


0.0041 (2.3) 


0.0021 (2.1) 


6 


0.00012 (1.0) 


0.00078 (0.74) 


0.00096 (2.1) 


0.00047 (2.2) 


7 


5.6E-05 (1.1) 


0.00048 (0.70) 


0.00022 (2.1) 


0.00015 (1.6) 



Table 5.9 
Convergence orders of estimated quantities for the convection-diffusion Cauchy problem ap- 
proximated using piecewise quadratic elements 7 = 0.001, 7;,^ = 10 (case 1) 
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N 


\\u — UfiW 


lk/.ll 


\uh\Sp + \Zh\Sa 


IIV(m-u^) -nW^if^Qf^ 


3 


0.020 (-) 


0.0014 (-) 


0.074 (-) 


0.12 (-) 


4 


0.034 (-) 


0.00028 (2.3) 


0.013 (2.5) 


0.11 (0.12) 


5 


0.026 (0.39) 


0.00011 (1.4) 


0.0025 (2.4) 


0.065 (0.76) 


6 


0.024 (0.12) 


8.3E-05 (0.4) 


0.00046 (2.4) 


0.043 (0.60) 


7 


0.023 (0.06) 


3.6E-05 (1.2) 


8.7E-05 (2.4) 


0.029 (0.57) 



Table 5.10 

Convergence orders of estimated quantities for the convection-diffusion Cauchy problem ap- 
proximated using piecewise quadratic elements 7 = 0.001, ■yi,^ = 10 (case 2) 



has some other interesting features. In particular for piecewise affine approxima- 
tion spaces the discrete solution can be shown to exist under very mild assumptions. 
Both symmetric stabilisation methods and the Galerkin least squares methods are 
considered in the analysis. Convergence of the method is obtained formally under 
abstract assumptions on the bilinear form that are shown to hold for three non-trivial 
examples. The actual performance of the method in practice depends crucially on 
the stability properties of the underlying PDE and when these are unknown, must 
be investigated numerically. Sometimes observed convergence orders are unlikely to 
match those predicted in Theorem 12.21 (except possibly for very small h), due to 
huge stability constants in the bound (|2.3p (cf. the Helmholtz equation for large wave 
numbers), or poor stability of the dual problem (cf. the Cauchy problem for Pois- 
son's equation). Another problem that may arise when ill-conditioned problems are 
considered, is poor conditioning of the system matrix. Even in the case of piecewise 
affine approximation the stabilisation corresponds to a very weak norm and in case 
the underlying problem is ill-posed this must be expected to show in the condition 
number. Clearly preconditioners for the linear systems arising is an important open 
problem. Other subjects for future work concerns the inclusion of hyperbolic prob- 
lems in the framework (see [5]) and the application of the method to data assimilation 
and boundary control. 
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